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Equation of state of a Fermi gas in the BEC-BCS crossover: a quantum Monte Carlo 

study 
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We calculate the equation of state of a two-component Fermi gas with attractive short-range in- 
terspecies interactions using the fixed-node diffusion Monte Carlo method. The interaction strength 
is varied over a wide range by tuning the value a of the s-wave scattering length of the two-body 
potential. For a > and a smaller than the inverse Fermi wavevector our results show a molecular 
regime with repulsive interactions well described by the dimer-dimer scattering length am = 0.6a. 
The pair correlation functions of parallel and opposite spins are also discussed as a function of the 
interaction strength. 

PACS numbers: 



Recent experiments on two-component ultracold 
atomic Fermi gases near a Feshbacli resonance have 
opened the possibility of investigating the crossover 
from a Bose-Einstein condensate (BEG) to a Bardeen- 
Gooper-Schrieffer (BGS) superfluid. In these systems the 
strength of the interaction can be varied over a very wide 
range by magnetically tuning the two-body scattering 
amplitude. For positive values of the s-wave scattering 
length a, atoms with different spins are observed to pair 
into bound molecules which, at low enough temperature, 
form a Bose condensate [l| . The molecular BEG state is 
adiabatically converted into an ultracold Fermi gas with 
a < and fci?|a| <C 1 PI, where standard BGS theory 
is expected to apply. In the crossover region the value 
of |a| can be orders of magnitude larger than the inverse 
Fermi wave vector kp and one enters a new strongly- 
correlated regime known as unitary limit @) Q • I'^ dilute 
systems, for which the effective range of the interaction 
Rq is much smaller than the mean interparticle distance, 
kpRn «C 1, the unitary regime is believed to be univer- 
sal |j,|a]. In this regime, the only relevant energy scale 
should be given by the energy of the noninteracting Fermi 
gas, 
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3 h^kj, 
10 m 
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The unitary regime presents a challenge for many-body 
theoretical approaches because there is not any obvi- 
ous small parameter to construct a well-posed theory. 
The first theoretical studies of the BEG-BGS crossover 
at zero temperature are based on the mean-field BGS 
equations |y]. More sophisticated approaches take into 
account the effects of fluctuations Q , or include explicitly 
the bosonic molecular field Q. These theories provide a 
correct description in the deep BGS regime, but arc only 
qualitatively correct in the unitary limit and in the BEG 
region. In particular, in the BEG regime the dimer-dimer 
scattering length has been calculated exactly from the so- 
lution of the four-body problem, yielding Um = 0.6a [9j. 



Available results for the equation of state in this regime 
do not describe correctly the repulsive molecule-molecule 
interactions 10). 

Quantum Monte Garlo techniques are the best suited 
tools for treating strongly-correlated systems. These 
methods have already been applied to ultracold degen- 
erate Fermi gases in a recent work by Garlson et al. [Il| ■ 
In this study the energy per particle of a dilute Fermi 
gas in the unitary limit is calculated with the fixed-node 
Green's function Monte Garlo method (FN-GFMG) giv- 
ing the result E/N = ^epo with £_ = 0.44(1). In a subse- 
quent work [ij , the same authors have extended the FN- 
GFMG calculations to investigate the equation of state 
in the BGS and BEG regimes. Their results in the BEG 
limit are compatible with a repulsive molecular gas, but 
the equation of state has not been extracted with enough 
precision. 

In the present Letter, we report results for the equa- 
tion of state of a Fermi gas in the BEG-BGS crossover 
region using the fixed-node diffusion Monte Garlo method 
(FN-DMG). The interaction strength is varied over a 
very broad range from —6 < —1/kpa < 6, including 
the unitary limit and the deep BEG and BGS regimes. 
In the unitary and in the BGS limit we find agreement, 
respectively, with the results of Rcf. 11] and with the 
known perturbation expansion of a weakly attractive 
Fermi gas [l3|- In the BEG regime, we find a gas of 
molecules whose repulsive interactions are well described 
by the dimer-dimer scattering length a™ = 0.6a. Results 
for the pair correlation functions of parallel and antipar- 
allel spins are reported in the various regimes. In the 
BEG regime we find agreement with the pair correlation 
function of composite bosons calculated using the Bogoli- 
ubov approximation. 

The homogeneous two-component Fermi gas is de- 
scribed by the Hamiltonian 



H 




(2) 



where in denotes the mass of the particles, i,j,... and 
label, respectively, spin-up and spin-down parti- 



N 



N/2, N being the total number of 



cles and A^- 

atoms. We model the interspecies interatomic interac- 
tions using an attractive square-well potential: V{r) — 
— Vb for r < Ro, and V{r) ~ otherwise. In order to en- 
sure that the mean interparticle distance is much larger 
than the range of the potential we use uRq = 10^^, where 
n = fc|./(37r^) is the gas number density. By varying the 
depth Vq of the potential one can change the value of the 
s-wave scattering length, which for this potential is given 
by a = Ro[l-ta.n{KoRo) / (KoRo)], where K^ = mVo/h'^. 
We vary Kq is the range: < Kq < tt/Rq- For KoRo < 
tt/2 the potential does not support a two-body bound 
state and a < 0. For KqRq > 7r/2, instead, the scattering 
length is positive, a > 0, and a molecular state appears 
whose binding energy et is determined by the trascen- 

dental equation 

where K'^ = Kq 

corresponds to the unitary limit where \a\ = oo and 

eb = 0. 

In a FN-DMC simulation the function /(R, r) = 
■!/'t(R-)^(R, t), where ^(R, t) denotes the wave function 
of the system and ■!/'t(R-) is a trial function used for im- 
portance sampling, is evolved in imaginary time r — it/h 
according to the Schrodinger equation 



eb\m/h^Ratan{KRo)/{KRa) = -1, 
\£b\m/h^. The value Kq = tt/{2Rq) 



dfi-R,T) 



- i?V2,/(R,r) + i^VR[F(R)/(R,r)] 

+ [EL{R)-Eref]fiIl,T). (3) 

In the above equation R = (ri, ...,rjv), -El(R) = 
'0t(R-)~^^'0t(R-) denotes the local energy, F(R) = 
2-0t(R)~^Vr'0t(R-) is the quantum drift force, D — 
h /(2m) plays the role of an effective diffusion constant, 
and Eref is a reference energy introduced to stabilize 
the numerics. The energy and other observables of the 
state of the system are calculated from averages over the 
asymtpotic distribution function /(R, r -^ oo). To en- 
sure positive definiteness of the probability distribution 
/ for fermions, the nodal structure of ipT is imposed as 
a constraint during the calculation. It can be proved 
that, due to this nodal constraint, the calculated energy 
is an upper bound to the eigenenergy for a given sym- 
metry [1,5]. In particular, if the nodal surface of i/'t were 
exact, the fixed-node energy would also be exact. 

In the present study we make use of the following trial 
wave functions. A BCS wave function 

^BCSCR) ^ A{(j){riy)(j){r22')--(l){rN^Ni)) , (4) 

and a Jastrow-Slater (JS) wave function 



i>js(R-) = l[^iru') 



aU^ 



aU' 



(5) 

where A is the antisymmetrizer operator ensuring the 
correct antisymmetric properties under particle ex- 
change. In the JS wave function, Eq. (O, the plane wave 
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FIG. 1: Energy per particle in the BEC-BCS crossover. Solid 
symbols refer to results obtained with the trial wave function 
"4^303, open symbols refer to the ones obtained with ^pjs- The 
dot-dashed line is the expansion Q holding in the BCS region 
and the dotted line corresponds to the binding energy 6^/2. 
Inset: finite size effects in the unitary limit —l/kpa — 0. 



orbitals have wave vectors ka = 2TT/L{£axX+£ayy+f-azz), 
where L is the size of the periodic cubic box fixed by 
nL^ — N, and £ are integer numbers. The correlation 
functions 0(r) and ip{r) in Eqs. Q-® are constructed 
from solutions of the two-body Schrodinger equation with 
the square-well potential V{r). In particular, in the re- 
gion a > we take for the function (/)(r) the bound-state 
solution (t>bsii") with energy e^ and in the region a < 
the unbound-state solution corresponding to zero scat- 
tering energy: (f>us{r) = {Rq - o) sin(i4:or)/[r sin(i4:oi?o)] 
for r < Rq and 4'usi'>') == 1 ^ o,/t for r > Rq. In the 
unitary limit, \a\ ^ oo, (j)bsir) ^ 0ns (r). 

The JS wave function tpjs, Eq. (O, is used only in 
the region of negative scattering length, a < 0, with a 
Jastrow factor ip{r) = 4'us{'f) for r < R. In order to 
reduce possible size effects due to the long range tail of 
(f'usir), we have used (p{r) — C1 + C2 exp(— ar) for r > R, 
with R < L/2 a matching point. The coefficients Ci and 
C2 are fixed by the continuity condition for (p{r) and its 
first derivative at r = .R, whereas the parameter a > 
is chosen in such a way that if{r) goes rapidly to a con- 
stant. Residual size effects have been finally determined 
carrying out calculations with an increasing number of 
particles N = 14, 38, and 66. In the inset of Fig. ^ we 
show the dependence of the energy per particle E/N on 
N in the unitary limit. Similar studies carried out in the 
BEC and BCS regime show that the value A^ = 66 is 
optimal since finite-size corrections in the energy are be- 
low the reported statistical error in the whole BEC-BCS 
crossover. We have also checked that effects due to the 
finite range Rq of the potential are negligible. 

The FN-DMC energies for TV = 66 atoms and the 
potential V{r) with uRq = lO^'' are shown in Fig. Q] 
and in Table ^ as a function of the interaction parame- 
ter — 1/kpa. The numerical simulations are carried out 
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FIG. 2: Energy per particle in the BEC-BCS crossover with 
the binding energy subtracted from E/N . SoUd symbols: re- 
sults with %I)BCS, open symbols: results with ^pJS- The dot- 
dashed line is as in Fig. and the dashed line corresponds to 
the expansion Q holding in the BEG regime. Inset: enlarged 
view of the BEG regime —l/kpa < —1. The solid line corre- 
sponds to the mean-field energy [first term in the expansion 
0], the dashed line includes the beyond mean-field correction 
[Eq. ©]. 



TABLE I: Energy per particle and binding energy in the BEG- 
BGS crossover (energies are in units of efg). 

-l/kpa E/N efc/2 E/N~eb/2 



-6 


-73.170(2) 


-73.1804 


0.010(2) 


-4 


-30.336(2) 


-30.3486 


0.013(2) 


-2 


-7.071(2) 


-7.1018 


0.031(2) 


-1 


-1.649(3) 


-1.7196 


0.071(3) 


-0.4 


-0.087(6) 


-0.2700 


0.183(6) 


-0.2 


0.223(1) 


-0.0671 


0.29(1) 
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0.42(1) 


0.2 


0.62(3) 
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0.87(1) 
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0.94(1) 
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both with the BCS wave function, Eq. Q), and with the 
JS wave function, Eq. 0. For —1/kpa > 0.4 we find 
that -ipjs gives lower energies, whereas for smaller val- 
ues of —1/kpa, including the unitary hniit and the BEC 
region, the function ipBCS is preferable. This behavior 
reflects the level of accuracy of the variational ansatz for 
the nodal structure of the trial wave function. We be- 
lieve that in the intersection region, —1/kpa ^ 0.4, both 
wave functions ipscs and ipjs give a poorer description 
of the exact nodal structure of the state, resulting in a 
less accurate estimate of the energy. In the BCS region, 
— 1/kpa > 1, our results for E/N are in agreement with 



FIG. 3: Pair correlation function of parallel, g2 (r), and (in- 
set) of antiparallel spins, gl {''')y for —l/kpa — (unitary 
limit), -l/kpa = -4 (BEG regime), -l/kpa = 4 (BGS 
regime) and for a noninteracting Fermi gas (FG). The dot- 
dashed line corresponds to the pair correlation function of a 
Bose gas with am = 0.6a and -l/kpa = —4 calculated using 
the Bogoliubov approximation. 
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In the unitary limit we find E/N = ^epc, with ^ = 
0.42(1). This result is compatible with the findings of 
Refs. [ll|, uM obtained using a different trial wave func- 
tion which includes both Jastrow and BCS correlations. 
The value of the parameter /3 = ^ — 1 has been measured 
in experiments with trapped Fermi gases [^ |3| , but the 
precision is too low to make stringent comparisons with 
theoretical predictions. In the region of positive scatter- 
ing length E/N decreases by decreasing kpa. At approx- 
imately — 1/kpa ~ —0.3, the energy becomes negative, 
and by further decreasing kpa it rapidly approaches the 
binding energy per particle eb/2 indicating the formation 
of bound molecules 13]. The results with the binding 
energy subtracted from E/N are shown in Fig. |3 In the 
BEC region, -1/kpa < -1, we find that the FN-DMC 
energies agree with the equation of state of a repulsive 
gas of molecules 



E/N - eb/2 



<^FG 
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IStt 



--{kpa„ 



^3/2, 



(7) 

where the first term corresponds to the mean-field en- 
ergy of a gas of molecules of mass 2m and density n/2 
interacting with the positive molecule-molecule scatter- 
ing length am, and the second term corresponds to the 
first beyond mean- field correction [la| . If for am we use 
the value calculated by Petrov et al. Q Om = 0.6a, we ob- 
tain the curves shown in Fig. |21 If, instead, we use am as 
a fitting parameter to our FN-DMC results in the region 



— l/kpa < —1, we obtain the value am/a = 0.62(1). A 
detailed knowledge of the equation of state of the homo- 
geneous system is important for the determination of the 
frequencies of collective modes in trapped systems 17], 
which have been recently measured in the BEC-BCS 
crossover regime [l8j . 

In Fig. 131 we show the results for the pair correlation 
function of parallel, §2 (r), and antiparallcl spins, 52 (^)- 
For parallel spins, §2 (r) must vanish at short distances 
due to the Pauli principle. In the BCS regime the ef- 
fect of pairing is negligible and §2 (r) coincides with 
the prediction of a noninteracting Fermi gas g] (r) — 
1 — 9 / {kpr)'^[sm{kpr) / kpr — cos{kpr)]'^ . This result con- 
tinues to hold in the unitary limit and it shows that in 
this regime the Fermi wave vector is not renormalized due 
to interactions, in agreement with the theory of Landau 
Fermi liquids jlQj . In the EEC regime the static structure 
factor S{k) of composite bosons can be estimated using 
the Bogoliubov result: S{k) = h'^k'^ /[2Muj{k)], where 
uj{k) = {h'^k'^ /AM'^ + gum^^k"^ /My/'^ is the Bogoliubov 
dispersion relation for particles with mass M = 2m, den- 
sity rim = n/2 and coupling constant g = AirTi Um/M. 
The pair distribution function 52 (f) of composite bosons, 
obtained through g2{r) ^ 1 + 2/NJ2k[S{k) - ^e-'^" 
using the value am = 0.6a, is shown in Fig. 13 for 

— l/kpci — —4 and compared with the FN-DMC re- 
sult. For large distances r ^ am, where Bogoliubov ap- 
proximation is expected to hold, we find a remarkable 
agreement. This result is consistent with the equation of 
state in the BEC regime and shows that structural prop- 
erties of the ground state of composite bosons are de- 



scribed correctly in our approach. For antiparallcl spins, 
gl (r) exhibits a large peak at short distances due to 
the attractive interaction. In the BEC regime the short 
range behavior is well described by the exponential decay 
gl {r) oc exjp{—2ry^]eb\m/h)/r^ fixed by the molecular 
wavefunction (j)ba{r). In the unitary regime correlations 
extend over a considerably larger range compared to the 
tightly bound BEC regime. In the BCS regime the range 
of <72 ('") i^ rnuch larger than kp and is determined by 
the coherence length S^q = h kp/imA), where A is the 
gap parameter. In this regime the wavefunction we use 
does not account for pairing and is inadequate to inves- 
tigate the behavior of gl (r). 

In conclusion, we have carried out a detailed study of 
the equation of state of a Fermi gas in the BEC-BCS 
crossover using FN-DMC techniques. In the BCS regime 
and in the unitary limit our results are in agreement with 
known perturbation expansions and with previous FN- 
GFMC calculations jll [13, respectively. In the BEC 
regime, we recover the equation of state of a gas of com- 
posite bosons with repulsive effective interactions which 
are well described by the molecule-molecule scattering 
length am = 0.6a recently calculated in Ref. jP|. 
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